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ABSTRACT 


Wireless sensors can be worn on soldiers or installed on vehicles to form 
distributed sensor networks to locate the source of sniper fire. A two-step source 
localization process is proposed for this sniper detection task. The time difference of 
arrival (TDOA) for the acoustic signals received by the sensors is first estimated using 
the generalized cross correlation (GCC) method. The estimated TDOA values are then 
used by the hybrid spherical interpolation/maximum likelihood (SI/ML) estimation 
method to estimate the shooter location. A simulation model has been developed in 
MATLAB to study the performance of the hybrid SI/ML estimation method. A wireless 
sensor network is simulated in NS-2 to study the network throughput, delay and jitter. 
Simulation results indicate that the estimation accuracy can be increased by increasing 
the number of sensors or the inter-sensor spacing. The constraint of small inter-sensor 
spacing on wearable sensors is found to degrade the estimation accuracy, but vehicular 
configuration providing larger inter-sensor spacing can help improve the estimation 
accuracy. The sensor topology should be well represented in all three dimensions to 
obtain desired estimation accuracy. The estimation accuracy is not adversely affected by 
sensor node failures or location perturbations. The NS-2 simulation results indicate that 
the wireless sensor network has low delay and can support fast information exchange 


needed in counter-sniper applications. 
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EXECUTIVE SUMMARY 


In urban warfare, sensor networks can be effectively used to detect, track, and 
monitor targets of interest. Sniper localization through the collaborative effort of multiple 
wireless acoustic sensors is of interest in this study. Wireless sensors worn on soldiers or 
installed on vehicles can be deployed to form a distributed sensor array to locate the 
source of sniper fire. A two-step source localization process is proposed in this work for 
counter-sniper applications: the time difference of arrival (TDOA) values are first 
determined using the generalized cross correlation (GCC) method and then the TDOA 
values are used by a hybrid spherical interpolation/maximum likelihood (SI/ML) 
estimation method to determine the shooter location. A simulation model has been 
developed in MATLAB to study the performance of the hybrid SI/ML estimation 
method. A wireless sensor network is simulated in NS-2 to study the network throughput, 


delay and jitter of IEEE WPAN and WLAN networks. 


The effects of number of sensor nodes, inter-sensor spacing, topology, noise, 
source-sensor distance, sensor node failure, location perturbation on the location 
estimation performance have been studied through MATLAB simulation. Simulation 
results show that the error variance can be reduced by increasing the number of sensor 
nodes or the inter-sensor spacing. To increase the location estimation accuracy, the 


number of sensors as well as the inter-sensor spacing should be increased. 


For the case of wearable sensors, the constraint of small inter-sensor spacing on 
the body degrades the estimation performance. On the other hand, vehicular 
configuration providing larger vertical and horizontal inter-sensor spacing can help 
improve the location estimation accuracy. The sensor node topology should be well 


represented in all three dimensions to obtain desired location estimation performance. 


The source location estimation accuracy is not adversely affected by sensor node 
failures or location perturbations. Results show that sensor node failures lead to a gradual 
degradation of accuracy and not a sudden total system failure. This demonstrates the 


advantage of using distributed wireless sensor networks to provide greater system 


XVil 


robustness. It is also observed that location perturbations do not adversely affect the 
location estimation. This illustrates the flexibility of the network configuration and 


demonstrates the suitability for potential field deployment. 


Network simulation results based on the NS-2 package show that the wireless 
sensor network can support a throughput between 20 kbps and 100 kbps, a low average 
delay of less than 10 ms and a low average jitter in the range of 5 ms to 10 ms. These are 
suitable for fast response network-based sniper detection applications that require low 


delay and fast data exchange. 
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I. INTRODUCTION 


Low-cost, low-power, multi-functional wireless sensor networks are becoming 
readily available. The sensor nodes are small in size and are able to sense, process data, 


and communicate with each other wirelessly. 


Wireless sensor networks have several potential applications, especially in the 
area of military operations [1][2]. In urban warfare, sensors can be deployed in buildings 
that have been cleared to prevent reoccupation. Sensors can be used to detect physical 
phenomena, such as motion, acceleration, acoustic and light, and the collected and 
processed sensor data can be transmitted through wireless links and displayed on mobile 


devices carried by soldiers. 


In urban warfare, sensor networks can be effectively used to detect, track, and 
monitor targets of interest. Sniper localization through the collaborative effort of multiple 
wireless acoustic sensors is of interest in this study. An array of acoustic sensors, 
microphones or piezoelectric pressure sensors, are deployed to detect and process 
supersonic shockwaves and muzzle blast produced by the gun shot [3]. The measured 
acoustic data are further processed using source localization techniques in order to 


estimate the location of the shooter. 


A. THESIS OBJECTIVE 

This thesis focused on source localization using wireless sensor networks. 
Wireless sensors worn on soldiers or installed on vehicles can be deployed to form a 
distributed sensor array to locate the source of sniper fire. To improve the accuracy of the 
source localization, a large number of sensors can be deployed to cover a large area of 
interest [4]. As illustrated in Figure 1, a distributed wireless sensor network, forming a 


large array of acoustic sensors, is used to detect and locate sniper fire. 


A two-step source localization process is proposed in this work for counter-sniper 
applications: the time difference of arrival (TDOA) values are first determined using the 
generalized cross correlation (GCC) method and then the TDOA values are used by a 


hybrid spherical interpolation/maximum likelihood (SI/ML) estimation method to 


1 


determine the shooter location. A simulation model is developed in MATLAB to study 
the performance of the hybrid SI/ML estimation method. A wireless sensor network is 
simulated in NS-2 to study the network throughput, delay and jitter of IEEE WPAN and 
WLAN networks. 


@ 
Sniper oes 
location Se 
OL O distributed 
aon wireless sensor 
bullet S, O 
Se Nig network 
(projectile) -.O O 
O-. OO 
O 
O O 
Figure 1. Sniper detection and localization using a distributed wireless sensor 
network. 


B. RELATED WORK 

Several physical phenomena can be used for sniper detection purposes. Pauli et al. 
have built the VIPER system that utilizes an infrared camera to detect the muzzle flash of 
the weapon [5]. It is augmented with a microphone to detect the muzzle blast for range 


estimation. Both sensors require direct line of sight. 


Commercial acoustic sniper detection systems, such as BBN’s Bullet Ears and 
Metravib’s PILAR, are also based on measurement of the time of arrival (TOA) of 
muzzle blasts and shockwaves produced by the gun shot. BBN’s Bullet Ears system 
utilizes one or two small arrays of microphones, providing estimates of range to the 


shooter and the caliber, speed and trajectory of the projectile [6]. 


Simon et al. have developed an experimental network-based counter-sniper 
system called PinPtr [4]. The system utilizes an ad hoc wireless sensor network built 
from inexpensive sensor nodes. After deployment, the sensor nodes synchronize their 


clocks, perform self-localization and wait for acoustic events. The sensors can detect 


muzzle blasts and acoustic shockwaves and measure their time of arrival (TOA). 
Utilizing a message routing service, the TOA measurements are delivered to a laptop 
computer, where the sensor signal processing algorithm calculates the shooter location 
estimate. The proposed solution utilizes time of arrival data of the measured shockwaves 
and muzzle blasts. From the measurements and the sensor positions a four-dimensional 
consistency function is formed. A quick search algorithm finds the maximum of this 


function. The location corresponding to the maximum is the shooter position estimate. 


Cc. THESIS OUTLINE 

The organization of the thesis is as follows: chapter II introduces the wireless 
sensor network as well as the various source localization techniques and sniper detection 
approaches. Chapter III describes the concept of source localization using wireless 
sensor networks. The two-step source localization technique proposed for counter-sniper 
applications will be discussed. Chapter IV presents the simulation model of the two-step 
source localization process, followed by the simulation results of the hybrid SI/ML 
estimation method using the MATLAB package simulating a sniper detection scenario. 
The throughput, delay and jitter of the wireless sensor network are simulated using the 
NS-2 network simulator. Chapter V concludes the thesis and highlights future work for 
further investigation. Appendix includes the MATLAB source codes used in the 


simulation studies. 
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Hl. INTRODUCTION TO WIRELESS SENSOR NETWORK, 
SOURCE LOCALIZATION AND SNIPER DETECTION 


This chapter discusses the various aspects of distributed wireless sensor networks, 
acoustic source localization and current sniper detection systems. We will introduce the 
concept of using wireless sensor networks to perform source localization for the purpose 


of sniper detection. 


A. DISTRIBUTED WIRELESS SENSOR NETWORKS 
Low-cost, low-power, multi-functional wireless sensor networks are becoming 
readily available. The wireless sensor nodes are small in size and are able to 


communicate with each other wirelessly by forming an ad hoc network. 


Wireless sensor networks are self-organizing networks with short-range 
communication and multi-hop routing. Their topologies frequently change due to 
mobility, fading and node failures. There are usually limitations on transmit power, 


memory, and computational capability in a sensor node [1]. 


The performance of a wireless sensor network generally depends on the following 
factors [1]: energy efficiency, latency, accuracy, fault tolerance, scalability, 
synchronization and localization, and throughput. Due to the interdependence of energy 
consumption, delay, and throughput, all these performance issues and metrics are tightly 


coupled. 


For this thesis research, the emphasis is not on any of the above performance 
metrics. Instead, we will focus on the issue of distributed array processing for 


determining the location of a signal source of interest. 


1. Architecture 

Wireless sensor networks can be organized in two architectures: flat or tiered. In a 
flat architecture, all sensor nodes are peers and are homogeneous in form and function. 
However, in a tiered architecture [7], as shown in Figure 2, sensor nodes form a hierarchy 
in which a sensor node at a given level performs a specific set of tasks on behalf of a 


subset of sensor nodes in the level below. 
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Figure 2. Tiered Architecture in Low Rate Wireless Personal Area Networks (LR- 
WPAN) (from [7]). 


In practice, very few sensor networks are entirely flat. Typically, a sensor network 
will connect to an internet via gateway nodes. In addition, sensor networks are often not 
physically homogeneous. For example, a network may become heterogeneous due to 
uneven battery drain across nodes. Tiered architectures are usually employed to take 
advantage of unevenly distributed resources by assigning resource-intensive roles to 


resource-rich nodes. 


2. Standards 

Efforts are under way to standardize the various layers of wireless sensor network 
communication protocols. IEEE 802.15.4 is the standard defined for Low Rate Wireless 
Personal Area Networks (LR-WPAN) [8], with the ZigBee Alliance [9] as its marketing 


and compliance certification organization. 


A LR-WPAN is a simple, low-cost communication network that allows wireless 
connectivity in applications with limited power and relaxed throughput requirements. The 
main objectives of a LR-WPAN are ease of installation, reliable data transfer, short-range 
operation, low cost, and a reasonable battery life, while maintaining a simple and flexible 


protocol suite. 


a. LR-WPAN Architecture 
The LR-WPAN architecture is specified in terms of a number of layers 
based on the open systems interconnection (OSI) seven-layer model, as shown in Figure 


3. Each layer is responsible for one part of the standard and offers services to higher 


Upper Layers 
$02.2 LLC 


layers. 






Figure 3. IEEE 802.15.4 LR-WPAN device architecture (from [8]). 


The IEEE 802.15.4 standard defines the physical layer (PHY) and the 
medium access layer (MAC). An IEEE 802.2 Type 1 logical link control (LLC) can 
access the MAC layer through the service specific convergence sub-layer (SSCS). The 
PHY layer specifies an operating frequency of 2.4 GHz with a basic bit rate of 250 kbps. 
There are alternate PHY specifications for 915 MHz and 868 MHz that operate at lower 
data rates of 40 kbps and 20 kbps, respectively. Details of the LR-WPAN modulation and 


spreading formats are summarized in Table 1. 









































PHY(MHz) | Frequency | Spreading parameters Data parameters 
band Chip rate | Modulation | Bit rate Symbol Symbols 
(MHz) (kchips/s) (kb/s) rate 
(ksym/s) 

868 868-868.6 300 BPSK 20 20 Binary 
915 902-928 600 BPSK 40 40 Binary 
2450 2400- 2000 O-QPSK 250 62.5 16-ary 

2483.5 Orthogonal 

Table 1. |. IEEE 802.15.4 LR-WPAN frequency bands and data rates (from [8]). 


b. LR-WPAN Network Topologies 


Depending on the application requirements, the LR-WPAN may operate in 


one of the two topologies shown in Figure 4: Star Topology or Peer-to-Peer Topology. In 


the Star Topology, the communication is established between devices and a single central 


controller, called the Personal Area Network (PAN) coordinator. Peer-to-Peer Topology 


also has a PAN coordinator; however, it differs from the star topology in that any device 


can communicate with any other device as long as they are within range of one another. 


Peer-to-Peer Topology allows more complex network formations to be 


implemented, such as mesh networking. A peer-to-peer network can be ad hoc, self- 


organizing and self-healing. It may also allow multiple hops to route messages from any 


device to any other device on the network. 
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Figure 4. Star and Peer-to-Peer topologies for a IEEE 802.15.4 LR-WPAN (from 
[8]). 


To form a larger network, supporting hundreds and thousands of sensor 
nodes, a single sensor cluster can reorganize and form a mesh of multiple neighboring 
clusters. Once predetermined application or network requirements are met, the PAN 
coordinator may instruct a device to become the cluster head of a new cluster adjacent to 
the first one. Other devices gradually connect and form a multi-cluster network structure, 
such as the one seen in Figure 5. The advantage of a multi-cluster structure is increased 


coverage area; the disadvantage is an increase in message latency. 





Figure 5. TEEE 802.15.4 LR-WPAN cluster tree network (from [8]). 
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B. ACOUSTIC SOURCE LOCALIZATION 

Sound localization can be accomplished by using observed differences in the 
sound signals received at different observation points, such as angle of arrival (AOA), 
time of arrival (TOA), time difference of arrival (TDOA), received signal strength (RSS) 
or steered response power (SRP), to estimate the direction and actual location of the 
sound source. Several methods, usually based on SRP and TDOA, have been developed 


for estimating an acoustical source location. 


Algorithms, like Steered Response Power-Phase Transform (SRP-PHAT) [10], 
are robust in the presence of reverberation. However, SRP-PHAT can be quite complex 
requiring the calculation of a large number of test points in the region of possible source 
locations. The location is chosen to be the point that produces the highest steered 


response power. 


Alternatively, the problem can be implemented as a two-step localization process 
[11]. In a two-step localization approach, the time difference of arrival (TDOA) is first 
determined by performing cross-correlation on the received signals [12]. The TDOA 


values are then used to estimate the location. 


There are generally three classes of TDOA estimators [12], including the general 
cross-correlation (GCC) approach, the maximum likelihood (ML) approach, and the 
phase transform (PHAT) approach. All these approaches attempt to calculate the cross- 
correlation in an optimal or suboptimal manner, and then select the time index of the peak 


of the result to be the TDOA estimate. 


The two-step localization algorithms are quite fast, but they lack robustness [12]. 
Frequently, errors do occur in the time delay estimates due to reflections of the sound 
source, which are sometimes greater in energy than the direct signal. The direct path can 
be obstructed or attenuated because of source and microphone directivity. Erroneous time 


delay estimates cause large errors in location estimation. 


For the counter-sniper application studied in this thesis, we propose the use of the 
two-step localization process due to its speed and the limited computational capabilities 


of the sensors. 
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Cc. SNIPER DETECTION 


A gun shot produces several detectable phenomena, such as muzzle blast and 


shockwaves [3] [6]. 


Muzzle blast is the acoustic signature associated with the ejection of the bullet 
from the sniper’s rifle. The muzzle blast is a loud, characteristic noise originating from 
the end of the muzzle and propagating spherically away at the speed of sound. The 
muzzle blast can be detected with acoustic sensors at ranges from several hundred meters 


out to more than a kilometer. 


Typical rifles fire projectiles at supersonic velocities, thereby producing acoustic 
shocks along their trajectory. Shockwaves can be detected acoustically at ranges from 
hundreds of meters out to more than a kilometer and can be used to accurately determine 
projectile trajectories. The shock waveform is distinctive and cannot be produced by any 


other natural phenomenon. 


As illustrated in Figure 6, the muzzle blast and shockwaves can be detected and 


processed by distributed wireless sensor networks equipped with acoustic sensors. 


bullet 





distributed 
wireless sensor 
network 
Muzzle wave front 0 Oe 2 
wave O = 
O O 
Figure 6. Shock wave front and muzzle wave detection by a group of distributed 


wireless sensors (after [4]). 


a 


As illustrated in Figure 7, wireless sensors may be mounted on vehicles for easy 
and fast deployment. The wireless sensors can also be worn on soldiers or integrated into 


their helmets. 





Figure 7. Wireless network-based counter-sniper system consisting of vehicular and 
wearable sensors. 


We have introduced the basics of wireless sensor networks, described the concept 
of source localization for acoustic signal sources and proposed the concept of using 
wireless sensor networks to detect sniper fire. In the next chapter, we will discuss in 
greater detail the various source localization techniques. We will discuss the two-step 


source localization technique proposed for counter-sniper applications. 
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HI. SOURCE LOCALIZATION USING WIRELESS SENSOR 
NETWORKS 


This chapter presents the details of source localization using distributed sensor 
networks. We will discuss the two-step source localization technique proposed for 


counter-sniper applications. 


Source localization is the process of determining the spatial location of a source 
based on multiple observations of the emitted signal source. The primary goal of source 
localization is accuracy. The accuracy of location estimates depends on a number of 
factors, including quantity and quality of sensors used, sensor placement relative to each 
other and the signal sources to be analyzed, ambient noise and reverberation, and number 
of active sources and their spectral content. The performance of localization techniques 
generally improves with the number of sensors in the array, especially when adverse 


environmental conditions are present. 


A. TWO-STEP LOCALIZATION PROCESS 
Figure 8 shows the schematic diagram of the two-step source localization 
technique considered in this thesis. The generalized cross correlation (GCC) technique is 


used to estimate the time difference of arrival 7,, for a pair of received signals m,(t) and 


m,(t). 


The Z,, values are then converted into range difference of arrival d,, values using 


the relationship 


A 


d,, = TD 


where 0 =345 m/s is the speed of sound [12]. Figure 9 shows the range difference of 
arrival at two acoustic sensors. The di, values along with the knowledge of the sensor 


positions are then used to generate hyperbolic curves which are then intersected to obtain 


a source location estimate. 


IS 


Figure 8. Two-step source localization model using GCC and hybrid SI/ML 
methods. Step 1 performs the TDOA estimation using GCC; Step 2 performs the 
location estimation in two stages: estimate determined by SI method in stage | is 

used as initial values for ML method in stage 2. 
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Figure 9. Range difference of arrival at two acoustic sensors. 


Location estimation can be performed using either a closed-form least- 


squares (LS) [13] or a maximum likelihood (ML) method [14]. The LS method is fast, 


and a closed-form solution can be obtained in a single iteration. The spherical 


interpolation (SI) method [13] is a widely used LS implementation in location estimation 


problems. In contrast, the ML method requires iterative gradient-descent search, and a 


good initial guess close to the solution is required in order to avoid divergence. 


However, the ML estimator is asymptotically unbiased and provides a more accurate 


estimate at low noise levels. 
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In this work, we propose the use of a hybrid SI/ML algorithm [14] for 
location estimation using the two-stage process, as illustrated in Figure 8: First, we 


calculate the source location §,, using spherical interpolation (SI) method. The estimate 
of the source location §,, is then used as initial value for the calculation of source 
location $,,, using the maximum likelihood (ML) method. A hybrid SI/ML algorithm has 


been proposed by Ajdler et al. [14] to perform acoustic source localization in a small 
room of size 3m X5 m X 3 m. In this thesis, a similar hybrid SI/ML algorithm is applied 


for acoustic source localization in an outdoor environment instead to detect sniper fire. 


At low noise levels, the SI method provides a rough estimate quickly. The 
ML method is computationally demanding and requires the use of an iterative technique, 
such as Gauss-Newton and Levenberg-Marquardt [15]. The convergence speed of the ML 
method can be increased by using the rough estimate from the SI method as initial 
estimate. Additionally, in situations where the ML method fails to converge, especially at 


high noise levels, the SI method can still provide a reasonable estimate. 


Other closed-form estimation methods like spherical intersection (SX) 
[13], planar intersection (PX) [13], and clustering-based techniques [16] can also be used 
to produce the initial guess of the source location. However, the SI method has greater 
noise immunity than the SX method and the SI method has lower error variance but 
slightly higher bias than the PX method. Clustering-based techniques are computationally 
less demanding than the SI method but are only applicable to sources that are very far 


away. 


Alternatively, the ML optimization problem can be solved by using 
alternating projection method [17] to avoid exhaustive multi-dimensional search. 
However, this is computationally more intensive than the closed-form methods and there 


is still no guarantee that the convergence to a global minimum can be obtained. 


ie 


1. Time Difference of Arrival Estimation 


Let m,(t) and m,(t) represent the time-delayed signals of a source signal received 
at two sensors as shown in Figure 9. If M,(@) and M,(@)represent the Fourier 
transforms of received signals m,(t) and m,(t), respectively, the TDOA 7Zcan be 


estimated using the GCC given as [12] 


@=argmax |” W (@)M,(@)M;(ae!do 
where 7 is an estimate of the delay between original source signal and the two sensors, 


and W (@)is the GCC weighting function. 


Common GCC weighting functions include maximum likelihood (ML) and phase 


transform (PHAT) [12]: 











ee IM, (@)||M, (a) 
“~~ |v, (@)| |, (of +|N, (@)| |M, (a) 
1 
W war (a) = IM, (@)M:; (a) 


where N(@) and N2(q@) are the estimated noise spectral density for the first and second 


sensors, respectively. 


Phase transform (PHAT) uses only the phase information derived from their 
respective signal power spectral density [12]. This phase-only procedure has the 
whitening effect on the signals and emphasizes primarily the channel effect. The PHAT 
has a peak at the relative time delay of the two signals and is used to estimate the TDOA 


information. 


The ML weights require knowledge about the spectrum of the sensor-dependent 
noises. The PHAT weights do not require this knowledge, and hence has been employed 
more often due to its simplicity. We choose the GCC-PHAT method due to its ability to 


reduce the degradation due to reverberation [10]. 
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2; Least-squares Spherical Interpolation (SI) Method 

The Spherical Interpolation method [13] is chosen for the closed-form least- 
squares (LS) source location estimation based on the measured range difference of arrival 
(RDOA) values. Figure 10 shows the notations used in SI method for a single source 


and two sensors i and the reference sensor 1. 

Let N denote the number of sensors. The vector of the spatial coordinates for the 
i-th sensor is denoted byx,=[x, y, z; li and the estimated position of the source is 
denoted byS=[x,  y, z,] . Sensor | is arbitrary selected as the reference node and 


thus the position of sensor | is set as the origin of the coordinate system (i.e., x, =0) 


The distance between the source and sensor 7 is denoted by the Euclidean norm 





and the distance from the origin to the point x; is denoted by R, = |x]. Similarly, the 


. The measured 














estimated distance from the origin to the source is denoted by R =| 


RDOA between sensor ij and sensor | is denoted by 
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sphere with radius R, 









source 
ee es sensor i 
‘ 
7 Xi 
D,=R, 
sensor | x 
(reference) _ xe 
= x, = 0 (origin) 
Figure 10. Notations used in Spherical Interpolation (SI) method for a single source, 


sensor i and the reference sensor | (from [13]). 


Assuming that sensor | is selected as the reference node (i.e., x, =), it can be 


shown [13] that the source location estimate is given as 


g=4(s’PiwPts) SPM WPLS. 





T 
where the projection matrix P;; = — ra : 
R -_ ad, 
o= RS ~ 45 ; 
Ry -_ dy 
d», 
d,, 
the range difference of arrival vector,d=| <> |, 
dy 
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. : X3 3 Z 
the sensor location matrix, S =| . 7 


Xy Yn <n 


Wi Wo Win 
; : Wo, Wo 8 Woy 

and the weight matrix, W=| - ; : : 
Wyi Wwo Wyn 


The range estimate is given as the Euclidean norm from the origin at sensor | to 
the source 


SJ=ya tyr +z, 


and the bearing estimate cos Q. is defined as the direction cosines from the origin at 


sensor | to the source 





cosQ. = 














The bearing estimate cos Q, is equivalent to the unit vector of the source location 


estimate §, as shown in Figure 11. 

















. source 
R= § Cesare 
s 
(0, 0, 0) 0” cos Q (unit vector) 
Figure 11. 


Range and bearing estimate of source location. 
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3. Maximum Likelihood Estimator 


The errors 77 in the measurement of the range difference of arrival (RDOA) d, 


are assumed to be zero mean Gaussian and independent for each sensor. We define X as 
the covariance matrix of the noise vector, Aas the vector of all the exactly measured 


RDOAs, andI‘as the vector containing the noisy measurements: 
P=A+7. 


Assuming that sensor | is chosen as the reference node, it can be shown [14] that 


the likelihood function of I’ given the source location s is 


f(T, s)= tron a5 


CH ipp 








where the RDOA error vector h(s) is represented by 


|s-»|-[s-al-4. 


h(s)= Is-al-Is-all-de (1.2) 











| || a |e | 


the noise covariance matrix is represented by 


Rl_ 
Re le 
Nie Ne|- 


ae? a Cs (1.3) 
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and the ML cost function is the exponent 


Jn (svat (sjz als). (1.4) 


From Equation (1.1) and Equation (1.4), it can be seen that maximizing the 


likelihood function f (I's) is equivalent to minimizing the ML cost function J,,, (s). 


Hence, the source location estimate 5,,, is obtained by minimizing the ML cost function 
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Sy, = argmin Jy, (5). 
s 


The cost function J,,(s) can be minimized using standard numerical 


optimization methods. The Levenberg-Marquardt method [15] provided by MATLAB is 


used to solve the minimization problem. 


4. Cramer Rao Lower Bound 

The minimum mean-squared error for any estimate of a non-random parameter is 
given by the Cramer Rao Bound [18]. The Cramer Rao Bound (CRB) gives a lower 
bound on the error variance of any unbiased estimate. The variance of any unbiased 


estimator of sis bounded by 
E| (8-s)(8-s)' |2F"(s) (1.5) 


where F (5 ) is the Fischer information matrix, expressed as 


ri)ee [eet amare] 1.6 





where Tis the vector containing the noisy RDOA measurements and f(I,s) is the 


likelihood function of I’ given the source location s. 


From Equation (1.1) and Equation (1.6), it can be shown that [14] 











where J(s) is the Jacobian matrix 
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From Equation (1.5) and Equation (1.7), it can be shown that the minimum 
variance of the estimates of the source position is equal to the trace of the inverse of the 


Fischer information matrix [14] 
min E| ($—s)(S-s)' |=0r(F"(s)). (1.8) 


Each diagonal element of the inverse Fischer information matrix corresponds to 
the minimal variance of one coordinate of the source position. The total variance on the 
estimation of the source position is the trace of the inverse Fischer information matrix. It 
can be shown [14] that the bound on the variance of the source location estimate is 


directly proportional to the noise variance. 


This chapter discussed the two-step source localization technique proposed for 
our counter-sniper applications. The generalized cross correlation scheme used to 
determine the time difference of arrival was briefly discussed and the least-squares 
spherical interpolation (SI) and maximum likelihood (ML) techniques have been 
summarized. In the next chapter, we will present the simulation results of the hybrid 
SI/ML estimation method as well as study the network throughput, delay and jitters of the 


wireless networks. 
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IV. SIMULATION RESULTS AND ANALYSIS 


In this study, the two-step source localization process, as discussed in Chapter II, 
is used to estimate the sniper location through the use of wireless sensor nodes forming a 
distributed sensor network to detect sniper fire. The generalized cross correlation (GCC) 
technique, as shown in Figure 8, is first used to estimate the time difference of arrival 
(TDOA) of the received acoustic signals at the sensors. The TDOA values, as shown in 
Figure 9, are then converted into range difference of arrival (RDOA) values to be used in 


the hybrid SI/ML estimation method to estimate the source location. 


This chapter presents the simulation model of the two-step source localization 
process, followed by the simulation results of the hybrid SI/ML estimation method using 
the MATLAB package simulating a sniper detection scenario. The MATLAB code used 
for the simulation is provided in Appendix. The throughput, delay and jitter of the 


wireless sensor network are simulated using the NS-2 network simulator. 


A. SIMULATION MODEL 

A simulation model is developed in MATLAB to study the two-step localization 
process. Figure 12 shows the simulation block diagram for the study of the two-step 
localization process. The study of the TDOA estimation using the GCC technique is 
carried out by examining the received signals at the sensors. The study of the hybrid 
SI/ML estimation is carried out separately by representing TDOA values as RDOA and 
modeling RDOA estimation errors as Gaussian noise [14]. The acoustic source is not 


explicitly simulated in this work. 





! ' estimation Ms 
oe ! 'TDOA ti 
oe ;converted — 
! ito RDOAd,, ' Stage 1 Stage2 | 
hamden aasios ! Step 2 ! 
Figure 12. Simulation block diagram for the study of the two-step localization 


process (based on Figure 8 in Chapter III). 
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1. Time Difference of Arrival Estimation 

In the simulation model, based on the known positions of the source and sensors, 
time-delayed copies of a source signal are generated for each sensor. We now consider 
the scenario of two sensors (see Figure 9). Figure 13 shows the simulated received 
signals at sensors 1 and 2. The received signal at sensor 2 is the same as the signal 


received at sensor | but is time delayed by 7,,= 40 ms. 


The TDOA between the two sensors is determined using the GCC method as has 
been discussed in Chapter III]. By applying the GCC technique on the two received 


signals, the TDOA estimate 7,, is determined to be 40.6 ms, which is close to the actual 


value of 7,, =40 ms. 





TDOA ,,= 40.0 ms 
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Figure 13. Received signals at two sensors. The received signal at sensor 2 is same as 
the signal received at sensor | but is time delayed by 7,,= 40 ms. 


Ze Hybrid SI/ML Location Estimation 
To study the hybrid SI/ML location estimation method proposed in the two-step 


localization process for counter-sniper applications, the TDOA information 7, 
determined using the GCC technique is now converted into RDOA values d, using the 


relationship d, =7,,V where V=345 m/s is the speed of sound [12]. In the simulation 
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study, however, the GCC method is not used. Instead, we have chosen to generate the 


measured RDOA values using an additive white Gaussian model as described below. 


The true RDOA dj values can be determined using the known source and sensor 
positions that are defined in the simulation scenario. Figure 14 shows the range difference 


of arrival for a single source and two sensors; sensor | is chosen to be the reference 


sensor. The locations of the source, sensor 2, and sensor | are (x,,y,,Z,), (%. 25%) and 
(0,0,0), respectively. The distances D,and D, are the Euclidean norms between the 


source and sensor | and the source and sensor 2, respectively. The RDOA d,, between 


sensor | and sensor 2 is given by 
d,, =D,—D, =7,,0 


is the magnitude of the difference between the time of arrival 1, at 





where 7,, = It, fe 


sensor | and the time of arrival ¢, at sensor 2. 
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Figure 14. Range difference of arrival d,, for a single source, sensor 2 and reference 


sensor | located at (x,,y,,z,). (%5, 5%) and (0,0,0), respectively. 


The measured RDOA values d, are modeled as 


A 


d,=d,+n 


where n is a white Gaussian noise with zero mean and variance o- . Figure 15 illustrates 


the additive white Gaussian noise model used to generate the measured RDOA values 
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de. The true RDOA values d,,are determined using the known sensor and source 


positions as illustrated in Figure 14. 


d, =,V d, =d,+n 
“po 
n 
Figure 15. Modeling of range difference of arrival (RDOA). 


The CP values are then used in the hybrid SI/ML estimation method to estimate 


the source location. As discussed in Chapter III, the hybrid SI/ML estimation method 
consists of two stages. As shown in Figure 12, the first stage of the hybrid SI/ML method 
produces the location estimate using the SI method. The SI estimate is then used to 
provide initial values for the iterative ML method used in the second stage. In the 
simulation study, location estimation performance of both SI and ML are studied using 


the performance metrics described below. 


Monte Carlo simulation runs are carried out and the estimated source location 


coordinates §, are compared with the actual source location coordinates s to obtain the 


performance metrics of error mean 


error variance 


1< 
o=— _ 
>| .=s= ul) 
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and root mean square (RMS) error given by square root of the second moment 


2 2 
YaVH +o 
where M is the number of simulation runs. 


Figure 16 shows location estimates of 100 Monte Carlo simulation runs of the SI 
and ML methods in a scenario with the source at 100 m from the origin; the sensors are 
spaced 10 m apart. Table 2 shows the error mean £/, error variance o* and RMS error y 
for 24 and 72 sensor nodes for the SI and ML methods. The notations Xs_S/, Ys_S/, 
Zs_SI, Rs_SI and Bearing_SI refer to the individual X, Y, Z range and bearing readings, 
respectively, for the SI method used in the first stage of the hybrid SI/ML method. 
Similarly, the notations Xs_ML, Ys_ML, Zs_ML, Rs_ML and Bearing_ML refer to the 
individual X, Y, Z range and bearing readings, respectively, for the ML method used in 
the second stage. From Figure 16, the estimated locations do not coincide with the actual 
source location at (0,100,0). For the case of 24 sensor nodes, the error variance O° is 
large and the position estimates are widely scattered. For the case of 72 sensor nodes, 
however, the ML method produces estimation with a low RMS error of less than 3 m, and 


the position estimates are closer to the source at (0,100,0). 
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Figure 16. Simulation output of the location estimation using the SI and ML methods 


in a scenario with the source at 100 m from the origin and the sensors spaced 10 
m apart: (a) 24 sensor nodes and (b) 72 sensor nodes. The number of Monte Carlo 
runs is 100. 
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Performance Metrics Number of Nodes 
24 72 
Error Mean (Xs_S/ -1.021 -7.6879 
Ms, in meters |Ys_SI -2.0383 20.688 
Zs_SI -1.3407 -0.7423 
IRs_SI -1.8955 20.992 
Bearing_SI 2.8112 3.7754 
Errormean (|Xs_ML 0.007947 -0.76399 
My, in meters |Ys_ML -0.37713 1.6088 
Zs_ML -0.01764 -0.17491 
Rs_ML -0.32973 1.6182 
Bearing _ML 1.535 0.54189 
Error variance Xs_S/ 11.914 11.224 
Go. in square |Ys_SI 359.54 10.753 
meters IZs_SI 7.5425 4.1336 
IRs_SI 353.59 12.112 
'Bearing_SI 11.599 1.9837 
Error variance Xs_ML 5.4121 1.2701 
a. in square} Ys_ML 196.71 5.7233 
meters Zs_ML 2.828 0.11484 
Rs_ML 195.4 5.7962 
Bearing _ML 1.4399 0.29564 
RMS error ([Xs_S/ 3.5996 8.3862 
Ys, i meters |Ys_S/ 19.071 20.946 
Zs_SI 3.0561 2.1644 
IRs_SI 18.899 21.279 
Bearing_SI 4.4161 4.0296 
RMS error Xs_ ML 2.3264 1.3615 
Ym, 0 meters |\Ys_ML 14.03 2.003 
Zs_ML 1.6818 0.38136 
Rs_ML 13.983 2.9008 
Bearing ML 1.9484 0).76765 











Table 2. Mean “, variance 0”, and RMS value y of the position errors for the SI and ML 


methods in a scenario with the source at 100 m from the origin and the sensors 
spaced 10 m apart. The number of Monte Carlo runs is 100. 
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B. HYBRID SY/ML ESTIMATION METHOD SIMULATION 

The effects of the number of sensor nodes, inter-sensor spacing, topology, noise, 
source-sensor distance, sensor node failure, and location perturbation on the estimation 
performance of the hybrid implementation of the closed-form least-squares spherical 
interpolation (SI) method and the iterative maximum likelihood (ML) method are studied. 
The main objective of this portion of the simulation is to study suitable wearable and 
vehicular sensor node configurations for the counter-sniper application as well as to 
determine the required number of sensors and the inter-sensor spacing to accurately 


estimate a source of sniper fire located a few hundred meters away. 


Figure 17 shows the schematic of a sensor node topology using the three- 
dimensional XYZ Cartesian coordinate system. By setting different source and sensor 
positions, different simulation scenarios can be constructed. Note that for the sensor 
node topology with regular inter-sensor spacing, a simplified notation of KxXLXM_ is 


used to describe a topology that has K,LandM sensor nodes along the 


X,Y and Z plane, respectively. 


source Q 


(X,¥42%;) 





Figure 17. Schematic of the sensor topology using the three-dimensional XYZ 
Cartesian coordinate system. White and grey circles represent sensors that are 
placed in the two-dimensional XY and the three-dimensional XYZ planes, 
respectively. 
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i, Effect of Number of Sensor Nodes 

The effect of the number of sensor nodes on the location estimation performance 
is studied using the scenario of a source located 100 m from the origin; the sensors are 
spaced 10 m apart. Several topologies are formed using 9 (1.e., 3X13), 12 (1.e., 2X23), 
24 (1.e., 4X2xX3), 36 (.e., 6X2X3), 48 (Le., 8X23), 60 (i.e., 10X23) and 72 (ie., 
12x2x3) sensor nodes. A 3xX1xX3 topology represents a single row of three sensor 
columns with three sensors each; a 12xX2x3 topology represents two rows of twelve 


sensor columns with three sensors each. 


Figure 18 and Table 3 show the error mean //, error variance 0” and RMS error 
yfor different sensor node topologies. For the cases of less than 12 sensor nodes, the 


errors of the ML were very large and are not shown in the plots in Figure 18. 


It can be seen that as the number of sensor nodes is increased, the error variance 
o° decreases for both SI and ML methods used in the two stages of the hybrid SI/ML 
algorithm. For example, when the number of sensors increases from 24 to 48, the error 
variance o° for the Rs values in the SI stage decreases quickly from 354 m’ to 36 m 
(approximately). The error variance o” of the ML stage is lower than the first SI stage 
and decreases from 195 m? to 7 m* (approximately). It is observed that when the error 
variance o” of the SI method in the first stage is high (i.e., for the case of 9 and 12 sensor 


nodes), it led to highly inaccurate results for the ML method in the second stage. 


The effect of the bias introduced by the SI method is clearly seen when the 
number of sensors is very high (1.e., 48 sensors or more). However, the ML method in the 
second stage is tolerant to the bias introduced by the SI stage. For example, the error 


mean £/ for Rs values using the SI method increased from 9 m to 21 m (approximately) 
when the number of sensors increased from 48 to 72; however, the RMS error y for Rs 


values for the ML method remained below 3.1 m. 
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Figure 18. Effect of number of sensor nodes on the location estimation performance 
of the hybrid SI/ML method: (a) error mean £, (b) error variance o”,, and (c) 
RMS error vy. The readings for 9 and 12 sensors are not shown due to highly 
erroneous results. Results are based on 100 Monte Carlo runs. 
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of Nodes 9 12 24 36 48 60 72 
Error Xs_ ST |-0.99945  -0.04349 +1.021 —-1.9582 |-2.8292 |-4.7036 |-7.6879 
mean“, [Ys_SJ -26.868 |30.456 [+2.0383 5.0957 8.7197 {13.989 [20.688 
in meters ZS_S/ {1.4638 1.8704  |+1.3407  |+1.6174 -1.3639 -1.1141 |-0.7423 
Rs ST -23.144 ~=|-26.308 |+1.8955 [5.2157 8.8395 [14.155 [20.992 
\Bearing 
|_ ST 21.487 22.94 2.8112 [2.3622 (2.2773 |2.7098 3.7754 
Error Xs_ML_ -3.49E+08 1.28E+09|0.007947 |-0.77733}-0.52786|-0.5967_ |-0.76399 
Imean 4, Ys_ML_ |1.60E+08 5.49E+08/+0.37713 [2.9459 [1.5521 [1.4334 [1.6088 
in meters 2S_ML_ -1.91E+08 /7.18E+08}-0.01764 |-0.34894-0.18401/-0.1626 -+0.17491 
Rs ML 4.37E+08 |1.57E+09}+-0.32973 [2.9664 1.5599 {1.4398 |1.6182 
\Bearing 
| ML 10.497 23.621 {1.535 0.97759 |0.54588 (0.46911 0.54189 
Error Xs SJ 10.163 (7.3159 {11.914 12.586 |12.01 {11.316 {11.224 
variance |Ys_ SJ [2088.6 [2156.9 359.54 {75.292 34.693 {18.157 [10.753 
o., in Zs_ SI 22.946 23.819 {7.5425 (6.8466 5.3038 [4.4692 (4.1336 
square Rs_SI (1560.7 1591.6 (353.59 {75.748 35.639 {19.304 {12.112 
meters Bearing 
|_ ST 1970.4 2098.1 11.599 {1.8786 {1.8306 {1.9483 {1.9837 
Error Xs ML 1.52E+18 4.12E+16.4121 [2.6975 |1.0198 (0.80249 |1.2701 
variance |Ys ML (2.76E+17 |7.71E+17|196.71 {32.131 (6.962 [4.2166 5.7233 
o,, in (Zs_ML (4.79E+17 [7.32E+17/2.828 0.90724 (0.30799 0.15477 (0.11484 
square Rs ML 2.26E+18 |1.54E+18|195.4 32.201 (6.9973 (4.2598 [5.7962 
meters Bearing 
| ML 447.78 233.55 {1.4399 0.34273 0.20057 (0.17955 0.29564 
RMS error|Xs_ST 3.341 2.7051 [3.5996 [4.0523 4.4737 [5.7827  |8.3862 
7 it Ys_SI 53.014 — 55.538 19.071 {10.063 {10.523 {14.624 20.946 
meters Zs_ SJ 5.0089 5.2266 3.0561 (3.0761 2.6766 [2.3897 2.1644 
Rs SI 45.786 47.788 18.899 {10.146 10.667 [14.821 [21.279 
\Bearing 
Bey 49.316 (51.228 4.4161 [2.7311 2.6489 (3.0482 4.0296 
RMS error|Xs_ML_ |1.28E+09 2.12E+08|2.3264 [1.8171 {1.1395 {1.0764 {1.3615 
Vj; Wl Ys_ML_ 5.49E+08 9.20E+08|14.03 6.3883 3.0612 2.5043 2.883 
meters Zs_ML_ |7.18E+08 8.98E+08|1.6818 [1.0144 (0.58468 (0.42569 |0.38136 
Rs ML 1.57E+09 1.30E+09|13.983 (6.4032 |3.0709 [2.5165 {2.9008 
\Bearing 
| ML 29021 17.475 {1.9484 {1.1395 0.70609 (0.63215 0.76765 
Table 3. | Effect of number of sensor nodes on the location estimation performance of the 


hybrid SI/ML method: error mean //, error variance o*, and RMS error y. The 


number of Monte Carlo runs is 100. 
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Another scenario is simulated to represent a number of soldiers wearing two 
sensors spaced 1 m apart. The sensors on the soldier are assumed to be mounted one 
each on the foot and the helmet. Each soldier is stationed 10 m apart and the source is at 
100 m from the origin. This leads to a 3-D sensor node arrangement, and the scenario 
considered the sensor node arrangements in the range from 4xX2x2 to 12xX2x2. A 
4x2x2 indicates a total eight soldiers wearing two sensors each in two rows of four 


soldiers each. 


Figure 19 shows the error mean, error variance and RMS error for the different 
number of soldiers. It can be seen that the estimation accuracy in the Z direction is poorer 
than in the X and Y directions. For example, for the case of 5x22, the RMS errors for 
the Xs, Ys, and Zs readings at the ML stage are 3 m, 10 m and 14 m, respectively. The 
poorer accuracy in the Z direction is due to the constraint that the wearable sensors are 


only | m apart in the Z direction. 


It is observed that by deploying more sensor nodes, the accuracy of the estimation 
can be improved. For example, twelve soldiers with two sensors each (i.e., 6X22) can 
detect a sniper 100 m away with an RMS error y of 3 m, 8 m and 15 m for the Xs, Ys, 
and Zs readings, respectively, at the ML stage. However, when the number of soldiers is 
greater than eighteen, the larger bias introduced by the SI method in the first stage may 
lead to less accurate estimates for the ML stage, especially for the Zs values. For 
example, forty-eight soldiers (12x22) can detect a sniper 100 m away with an RMS 


error Y of 3 m, 4m and 17 m for the Xs, Ys, and Zs readings, respectively, using the ML 


method. 
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Figure 19. Performance of the hybrid SI/ML method for simulated wearable 


configuration (two sensors per soldier) with different number of soldiers: (a) error 
mean JZ, (b) error variance o°, and (c) RMS error 7. Results are based on 100 


Monte Carlo runs. 


Ze Effect of Inter-Sensor Spacing 


To study the effect of inter-sensor spacing on the location estimation 


performance, a scenario consisting of 24 sensor nodes (6X2X2) is simulated. The source 


is situated 100 m from the origin and the inter-sensor spacing is varied from 10 m to 50 


m. Figure 20 and Table 4 show the error mean //, error variance o” and RMS error y 


for different inter-sensor spacing. 


It is observed that when the inter-sensor spacing increases, the estimation 


: : 2 : : 
accuracy increases. The error variance o° decreases rapidly when the inter-sensor 


spacing is increased. For example, at the inter-sensor spacing of 10 m, the error variance 


o° of the Rs values obtained at the ML stage was 195 m’; at the inter-sensor spacing of 


20 m, the error variance o” quickly decreased to 7.6 m”. 
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Figure 20. Effect of inter-sensor spacing on location estimation performance of the 
hybrid SI/ML method : (a) error mean JZ, (b) error variance o”, and (c) RMS 


error Y. Results are based on 100 Monte Carlo runs. 
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Inter- 
sensor 
Spacing lm 5m 10m | 20m | 30m | 40m 50m 
IMean ([Xs_SI 1.4874 3.8611} -1.021) -3.9662) -5.0268} -6.0904| -7.0891 
error £l,, |\Ys_SI -96.052,  -56.341) -2.0383) 6.8829) 5.3421} 4.125) 2.4439 
in meters 259! 1.014 2.3276 -1.3407| -3.667| -5.3771| -7.7525| -10.657 
Rs_SI -92.607 -49.787| -1.8955| 7.0553} 5.6146 4.5986, 3.2445 
Bearing_ 
ST 75.399 42.603) 2.8112) 2.8987, 3.9867, 5.4073) 7.1224) 
Mean Xs_ML | -6.54E+07/4.01E+07| 0.00795] -1.0804| -0.6564|-0.30166| -0.04878 
error Ys_ML 1.11E+08)2.03E+08|-0.37713| 3.4316) 1.5148} 0.55075) 0.055739 
fy, 1 Zs_ML -2.82E+07 3.03E+06|-0.01764|-0.75852)-0.503 1 1|-0.24626| -0.04254 
meters ([Rs_ML 1.35E+08|2.46E+08|-0.32973| 3.4495) 1.5218} 0.55344) 0.056847 
Bearing_ 
ML 26.764 10.791) — 1.535] 0.91882) 0.57074} 0.34763) 0.2187 
Error Xs_SI 0.054935) 1.19E+01} 11.914; 5.0661) 2.1157) 0.99307 0.451 
variance |Ys_S/ 244.13) 1.97E+03) 359.54 14.286) 1.6202} 0.1256) 0.014497 
oy, in |Zs_SI 0.024743|5.78E+00) 7.5425} 3.6117) 1.5424, 0.7841) 0.41151 
square Rs_ST 208.03) 1.38E+03} 353.59, 15.188) 1.9432} 0.22131) 0.002205 
meters |Bearing_ 
ST 3060.8) 3297.6} 11.599, 1.9415} 0.96781} 0.4869} 0.25618 
Error = |Xs_ML 1.64E+16|7.94E+16, 5.4121) 1.2418) 0.44208) 0.21931) 0.11579 
variance |Ys_ML 4.77E+16,2.39E+17| 196.71} 7.6353} 0.97013} 0.20748} 0.032234 
o.,, in \Zs_ML 4.07E+15)1.57E+16, 2.828) 0.81316) 0.32197] 0.17749] 0.10389 
square IRs_ML 6.74E+163.16E+17| 195.4) 7.7062) 0.98182) 0.2098] 0.032316 
meters |Bearing_ 
ML 709.58 174.95) 1.4399) 0.27527] 0.12701} 0.055) 0.025056 
RMS  Xs_SI 1.64E+16/5.18E+00) 3.5996 4.5604) 5.233) 6.1714) 7.1209 
error Y,, |Y¥s_ST 4.77E+16,7.17E+01, 19.071) 7.8524) 5.4917) 4.1402) 2.4469 
in metersZ5_9/ 4.07E+15)3.35E+00, 3.0561) 4.1302) 5.5186) 7.8029 10.676 
Rs_ST 6.74E+166.21E+01} 18.899) 8.0601} 5.7851) 4.6226, 3.2449 
Bearing_ 
ST 709.58 71.503) 4.4161) 3.2162) 4.1063) 5.4522) 7.1404 
RMS Xs_ML 1.44E+08)/2.85E+08) 2.3264) 1.5521} 0.93432} 0.55705) 0.34377 
error Ys_ML 2. 45E+08)/5.29E+08} 14.03} 4.4058) 1.8069) 0.71471} 0.18799 
Vu, 0 Zs_ML 6.98E+07 1.25E+08} 1.6818) 1.1783) 0.75834} 0.48799} 0.32511 
meters [Rs_ML 2.93E+086.14E+08} 13.983) 4.4278) 1.816 0.7184) 0.18854 
Bearing_ 
ML 37.761 17.07, 1.9484, 1.0581} 0.67287] 0.41934, 0.26998 
Table 4. Effect of inter-sensor spacing on location estimation performance of the hybrid 


SI/ML method: error mean £/, error variance o’, and RMS error y. Results are 


based on 100 Monte Carlo runs. 
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Additional simulations are carried out to take into account the limited vertical 
inter-sensor spacing when the sensors are worn on soldiers or installed on vehicles. It is 
assumed that there is limited height of less than 2 m for the soldiers and less than 7 m for 


vehicles with elevated masts. 


a. Wearable Configuration 

To study the wearable configuration for the soldiers, another scenario is 
simulated to represent twelve soldiers (6X22) with two sensors per soldier for different 
Z-spacing values of 1.0 m, 1.25 m, 1.5 m, 1.75 m and 2 m. Figure 21 shows the RMS 


error ¥ for the different Z-spacing values. 


It can be seen that the RMS error decreases when the Z-spacing is 
increased. For example, the RMS errors for Zs readings at the ML stage decreased from 
15 m to 8 m when the Z-spacing was increased from 1 m to 2 m; however, the RMS 
errors for Xs and Ys readings decreased slightly from 3 m and 8 m to 2 m and 7 m, 
respectively. Hence, it is recommended that the two sensors worn on the soldiers should 
be separated as far apart as practically possible; however, considering the constraint of 
the average height of a person ranging from 1.5 m to 2 m, the assumption of 1 m for the 
Z-spacing separation, say between the head to the waist level, is a reasonable assumption 


for a wearable configuration. 
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Figure 21. RMS error for different values of Z-spacing for 6X22 sensor nodes, 


representing a wearable configuration. Results are obtained by averaging output 
from 100 Monte Carlo simulations. 
b. Vehicular Configuration 
A vehicular configuration can provide a larger vertical inter-sensor 
spacing compared to the wearable configuration. The second sensor can be mounted on 
an elevated mast to provide a larger sensor spacing in the Z direction and will improve 


the estimation accuracy. 


To study the vehicular configuration, the scenario representing eight 
vehicles (4X22) with two sensors installed vertically on each vehicle is simulated. The 
vehicles are spaced 10 m apart and the source is situated 100 m from the origin. The Z- 
spacing values are varied from 1 m to 7 m. Figure 22 shows the RMS error for different 


Z-spacing values. 


When the Z-spacing is increased, the estimation accuracy increases. For 


example, with a Z-spacing of 3 m, the RMS error y at the ML stage was 3 m, 14 m and 4 


m for Xs, Ys, and Zs readings, respectively; with a Z-spacing of 7 m, the RMS errors 


decreased to 2 m, 12 m and 2 m for Xs, Ys, and Zs readings, respectively. 
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Figure 22. RMS error y for different values of Z-spacing, from 1 m to 7 m, 


representing a vehicular configuration. Results are obtained by averaging output 
from 100 Monte Carlo simulations. 


A similar scenario consisting of eight vehicles is repeated, but the 
simulation is carried out with different values of inter-vehicle spacing, from 10 m to 40 


m. Figure 23 shows the RMS error for different vehicular spacing values. 


It can be seen that the estimation accuracy can be improved by increasing 
the inter-vehicle spacing beyond 10 m. For example, for an inter-vehicle spacing of 10 
m, the RMS error at the ML stage was 15 m for the Rs readings; for an inter-vehicle 


spacing of 40 m, the RMS error was significantly reduced to 4 m for the same reading. 


In this scenario, however, wireless links that can cover a longer physical 
range will be needed. For instance, the IEEE 802.11 WLAN network with an effective 
range of up to 100 m can be used instead of the IEEE 802.15.4 WPAN that has an 


effective range of up to only 10 m. 
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Figure 23. RMS error for different values of inter-vehicular spacing of 10 m, 20 m, 


30 m, and 40 m, representing a vehicular configuration. Results are obtained by 
averaging output from 100 Monte Carlo simulations. 


3. Effect of Topology 

After studying the effect of the number of sensor nodes and the inter-sensor 
spacing, further simulations are carried out to find out a suitable sensor topology for the 
counter-sniper applications and the required number of sensors to be worn on soldiers. To 
study the effect of topology on the location estimation performance, different two- 
dimensional (4xX2x1 and 8x2xX1) and three-dimensional topologies (4x2x2 and 
4x2x3) are simulated. Figure 24 shows four different simulated topologies used in the 
study. The 4x2x1 and 8x2x1 topologies refer to the layout of eight and sixteen sensors 
on the XY plane, respectively. The 4x2x2 and 4x2x3 topologies refer to the layout of 
eight sensor columns with two sensors and three sensors, respectively, along the Z axis. 


Each sensor is separated by a distance of 10 m. 
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Figure 24. Different sensor topologies: (a) 4X21, (b) 8X2x1, (c) 4X22 and (d) 
4x2x3. White and grey circles represent sensors that are placed in the two- 
dimensional XY and three-dimensional XYZ planes, respectively. 


Figure 25 shows the RMS errors for the four sensor topologies considered. 
Simulation results indicate that although the estimation errors are low in the X direction 


for the case of eight sensors lying on the XY plane (i.e., 4x21 topology), the estimation 
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errors in the Y and Z directions are relatively large. For example, RMS error of the Xs, Ys 


and Zs values at the ML stage for the 4x2x1 topology was 10 m, 40 m and 30 m, 
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Figure 25. RMS error y for different two-dimensional (42x 1 and 8x21) and 


three-dimensional topologies (4x22 and 4x2x3). Results are obtained by 
averaging output from 100 Monte Carlo simulations. 


The same observation of large estimation errors in the Y and Z directions can be 
seen when the number of sensors lying on the XY plane is increased from eight to sixteen 
(1.e., 8X2X1 topology). For example, the RMS error for the Xs, Ys and Zs values at the 
ML stage for the 8x21 topology was 15 m, 35 m and 30 m, respectively. However, it is 
observed that when the number of sensor nodes in the Z direction is increased (1.e., 
4x2x2 and 4x2x3 topologies), there is significant improvement in the estimation 
accuracy in the all three directions. For example, the 4x2x1 and 8x2x1 topologies 
produce the largest error of 30 m in the Zs readings at the ML stage while the 4x2x2 and 
4x2x3 topologies produce the smallest error of 3 m in the Zs readings. Hence, it is 
recommended that the sensor topology should be well represented in all three directions 


to obtain good estimation accuracy in all three dimensions. 


Another scenario is carried out to represent topologies for sensors worn on 


soldiers. For wearable sensors, it is assumed that there is a limited height of less than 2 m 
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where the sensors can be worn on the soldier’s foot and helmet. The simulation is carried 
out by having one topology of two wearable sensors spaced 1 m to 1.5 m apart and 


another topology of three wearable sensors spaced 0.5 m to 0.75 m apart. 


Figure 26 shows the RMS error yfor different number of vertical sensors per 


soldier: two wearable sensors spaced | m to 1.5 m apart (6X2X2) and three wearable 
sensors spaced 0.5 m to 0.75 m apart (6X2x3). It can be seen that the estimation 
accuracy is comparable for the two-sensor (6X2xX2) and three-sensor wearable 
configurations (6X23). For example, the topologies of 6X22 (Z-spacing = | m) and 
6xX2x3 (Z-spacing = 0.5 m) produced similar RMS errors of 3 m, 8 m and 15 m for Xs, 
Ys and Zs readings, respectively, at the ML stage. Similarly, the topologies of 6x22 (Z- 
spacing = 1.5 m) and 6X23 (Z-spacing = 0.75 m) produced similar RMS errors of 3 m, 
7 mand 11 m for Xs, Ys and Zs readings, respectively. Thus, it is recommended that two 


sensors per soldier will be sufficient. 
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Figure 26. RMS error yfor two-sensor (6X2X2 with Z-spacing of 1 m and 1.5 m) 


and three-sensor (6X23 with Z-spacing of 0.5 m and 0.75 m) topologies, 
representing wearable configuration. Results are obtained by averaging output 
from 100 Monte Carlo simulations. 


4. Effect of Noise 


Measurement noise is expected to have a significant effect on the estimation 


accuracy [14]. The variance in acoustic travel time o7 along a ray path from source to 


receiver is proportional to the propagation distance R, as given by [3] 


RL 
0, =0.4y° 
D 
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where Lrepresents the correlation length of atmospheric turbulence, v is the speed of 
sound and y’ is the mean-squared sound-speed fluctuations near the ground. However, 
the proportionality property is valid only for short distances of less than 60 m. A higher 
order approximation is required to provide a more realistic estimation. In this study, the 


RDOA noise variance is modeled as a second-order approximation 


where a is aconstant that is set to 0.01 for the study. 


To study the effect of noise on the location estimation performance, a scenario 
consisting of 24 sensor nodes (6X2X2), as shown in Figure 27, with a 10-m spacing 


between adjacent sensors is simulated. The source is placed 100 m from the origin. 


source 7) 


(0,100, 0) 





Figure 27. Simulated scenario of 24 sensor nodes (6X22) with a 10-m spacing 
between adjacent sensors is simulated. The source is placed 100 m from the 
origin. 


Figure 28 and Table 5 show the error mean //, error variance o°, and RMS error 


y for different noise levels. 
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Figure 28. Effect of noise on location estimation performance of the hybrid SI/ML 
method with different noise variance oO. for: (a) error mean //, (b) error variance 


o°, and (c) RMS error 7, representing a 6X22 topology. Results are based on 
100 Monte Carlo runs. 
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Noise 0.001 0.01 0.1 1 
variance oO. 
Error mean Xs_SI -2.5658 -2.5682 -2.5822 -1.021 
f,, inmeters | Ys_SI 11.133 11.141 11.115 -2.0383 
Zs_SI -2.1288 -2.1309 -2.1485 -1.3407 
Rs_SI 11.183 11.191 11.168 -1.8955 
Bearing_SI 1.7183 1.7198 1.7261 2.8112 
Error mean Xs_ML -2.1565 -1.4044 -1.4486 0.007947 
My, 1 meters | Ys_ML 10.818 9.0713 9.2698 -0.37713 
Zs_ML - 1.6793 -0.95559 -0.97037 -0.01764 
Rs_ML 10.854 9.0848 9.2844 -0.32973 
Bearing ML | 1.4118 0.89127 0.91088 L355 
Error variance | Xs_S/ 2.57E-05 0.002569 0.25342 11.914 
a. in square | _Ys_S/ 0.000386 0.038694 3.8292 359.54 
aieicee Zs_SI 1.84E-05 0.001845 0.18198 7.5425 
Rs_SI 0.000394 0.039445 3.9012 353.59 
Bearing_SI _ | 9.72E-06 0.000973 0.09733 11.599 
Error variance | Xs_ML 2.66E-01 0.04107 0.1096 5.4121 
ay in Ys_ML 1.837 1.1165 3.4264 196.71 
square meters Zs_ML 2.88E-01 0.02666 0.063778 2.828 
Rs_ML 1.8589 1.1228 3.4485 195.4 
Bearing_ML | 1.39E-01 0.017059 0.041712 1.4399 
RMS error Xs_SI 2.5658 2.5687 2.6308 3.5996 
y;, inmeters | Ys_SI 11.133 11.143 11.286 19.071 
Zs_SI 2.1288 2.1313 2.1905 3.0561 
Rs_SI 11.183 11.193 11.341 18.899 
Bearing_SI 1.7183 1.7201 1.7541 4.4161 
RMS error Xs_ML 2.2173 1.419 1.4859 2.3264 
Vy, meters | Ys_ML 10.902 9.1327 9.4528 14.03 
Zs_ML 1.7629 0.96944 1.0027 1.6818 
Rs_ML 10.939 9.1464 9.4683 13.983 
Bearing_ML | 1.4603 0.90079 0.9335 1.9484 





Table 5. 


Effect of noise on the location estimation performance of the hybrid SI/ML 





method with different noise variance o for: error mean /Z, error variance Oo”, 
and RMS error 7, representing a 6X22 topology. Results are obtained by 
averaging output from 100 Monte Carlo simulations. 


It can be seen that as the noise variance o* increases, the error variance o 


n 
increases. For example, the error variance o* for the Rs value at the ML stage increased 


from 1.1 m* to 3.4 m’ when the noise variance o> increased from 0.1 to 1. 


49 


Although it is expected that a noisier environment will lead to greater estimation 


errors, it is interesting to note that at high noise levels (i.e.,o7 =1), ML actually produces 
an error mean that is better than for low noise levels (i.e., Oo. = 0.001, 0.01 or 0.1). For 


example, the error mean for the Rs values for the ML stage are 11 m and -2 m for o- 


values of 0.1 and 1, respectively. This anomaly is likely due to the inaccurate SI value 
that is used as the initial value for the ML stage. At low noise levels, the SI stage 
introduces large bias of approximately 11 m and leads to similarly large errors of 
approximately 10 m for the ML stage. The smaller bias of 2 m introduced by SI at high 


noise levels actually leads to better estimation accuracy in the ML stage. 


a. Effect of Source-Sensor Distance 
To study the effect of the source-sensor distance on the location estimation 
performance, a scenario consisting of 24 sensor nodes (6X2X2) is simulated. The 


distance between the source and the sensor is varied from 10 m to 200 m. Figure 29 and 
Table 6 show the error mean, error variance 0°, and RMS error y for the various 


source-sensor distances in this range. 


As the source-sensor distance increases, the error variance o” increases. It is seen 
that the error variance of the SI method in the first stage increases much rapidly than that 
in the ML stage. For example, at the first SI stage, the error variance of Rs value 
increases from 4 m* to 354 m” when the source-sensor distance is increased from 50 m to 
100 m; the error variance at the ML stage increases from 2 m’ to 195 m’, 


correspondingly. 


However, when the source-sensor distance increases from 150 m and beyond, the 
ML method produces very divergent results due to the erroneous estimates provided by 


the SI method and are not shown in Figure 28. 
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Figure 29. 


(c) 


Effect of source-sensor distance on location estimation performance of the 


hybrid SI/ML method: (a) error mean JZ, (b) error variance o, and (c) RMS 
error Y, representing a6X2xX2 topology. Results are obtained by averaging 


output from 100 Monte Carlo simulations. Readings for 150 m and beyond are not 
shown due to erroneous results. 
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Source- 10m 50m 100 m 150m 200 m 
sensor 
distance 

Error mean Xs_SI 0.55436} -1.9831 -1.021 3.6574 7.7222 

f,, inmeters | Ys_SI -2.3537| 3.4415 -2.0383 -48.569 -112.68 
Zs_ST -1.6988] -1.8335 -1.3407 1.7466 4.6552 
Rs_SI -2.1474| 3.5276 -1.8955 -45.186 -99.574 
Bearing_ 13.165] 2.8987 2.8112 16.813 42.603 
SI 

Error mean Xs_ML -0.01412| -0.54358 | 0.007947 | 1.21E+08 | 2.76E+07 

My, in meters | Ys_ML 0.005905] 1.7556 | -0.37713 | 3.36E+08 | 4.72E+08 
Zs_ML -0.0134| -0.38831 | -0.01764 | -7.11E+07 | -1.56E+08 
Rs_ML 0.005981] 1.7648 | -0.32973 | 4.27E+08 | 5.73E+08 
Bearing_ 0.15755} 0.93126 1.535 5.4505 10.156 
ML 

Error variance | Xs_S/ 0.004201) 1.2665 11.914 36.192 47.508 

exe in square | Ys_SI 0.005741) 3.5715 359.54 3341.5 7866.5 

nieiere Zs_SI 0.001391) 0.90292 7.5425 18.249 23.135 
Rs_SI 0.005977] 3.7969 353.59 2707.9 5533.6 
Bearing_ | 0.000324) 1.9415 11.599 1240 3297.6 
SI 

Error variance | Xs_ML 0.000607} 0.32603 5.4121 | 2.18E+18 | 5.71E+17 

G,. in Ys_ML 0.000151 2.001 196.71 | 3.46E+18 | 1.32E+18 

square meters Zs_ML 0.000558] 0.21051 2.828 | 1.70E+17 | 1.96E+17 
Rs_ML 0.000149) 2.0198 195.4 | 5.75E+18 | 2.00E+18 
Bearing_ | 0.025758) 0.28883 1.4399 86.278 144.19 
ML 

RMS error Xs_SI 0.55813) 2.2802 3.5996 7.0404 10.351 

y,, inmeters | Ys_SI 2.3549) 3.9262 19.071 75.501 143.4 
Zs_SI 1.6992} 2.0651 3.0561 4.6151 6.6938 
Rs_SI 2.1487 4.03 18.899 68.918 124.29 
Bearing_ 13.165) 3.2162 4.4161 39.021 71.503 
rod 

RMS error Xs_ML 0.0284) 0.78836 2.3264 | 1.48E+09 | 7.56E+08 

Vy, in meters | Ys_ML 0.013616} 2.2546 14.03 | 1.89E+09 | 1.24E+09 
Zs_ML 0.027152] 0.60108 1.6818 | 4.18E+08 | 4.70E+08 
Rs_ML 0.013599) 2.2659 13.983 | 2.44F+09 | 1.53E+09 
Bearing_ 0.2249) 1.0752 1.9484 10.77 15.727 
ML 





Table 6. 


Effect of source-sensor distance on location estimation performance of the hybrid 


SI/ML method: error mean //, error variance o*, and RMS error 7, representin 
ra p g 


a 6X2x2 topology. Results are obtained by averaging output from 100 Monte 
Carlo simulations. 
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Another scenario is carried out to simulate the effect of source-sensor distance on 
the location estimation performance for a group of soldiers equipped with two sensors 
each. A scenario consisting of 24 sensors (6X2X2) is simulated to represent twelve 
soldiers who are individually equipped with two sensors that are spaced 1 m apart on 
their body. Each soldier is stationed 10 m apart, and the source is situated 50 m to 200 m 


from the origin. Figure 30 shows the RMS error y for the 6X2x2 soldier configuration 


with different source-sensor distances: 50 m, 100 m, 150 m and 200 m. 


It is seen that when the source-sensor distance increases, the RMS error vy 
increases. For example, at a source-sensor distance of 100 m, the RMS error yfor Xs, Ys 


and Zs readings at the ML stage is 3 m, 8 m and 15 m, respectively; at a source-sensor 
distance of 200 m, the RMS error yfor Xs, Ys and Zs readings at the ML stage is 8 m, 45 


m and 32 m, respectively. 




































































90 
80 
70 © 50m 
2 60 
© 50 w 100m 
= > 150m 
aa 20 oO 200m 
410 
0 
S S > yw y XY wy wy 
42” Ao7 A27 ee 7 Aue 42 of 12 7 AP? S no 7 Ra 
eg & 
Figure 30. RMS error yfor different values of source-sensor distances: 50 m, 100 m, 


150 m and 200 m, representing a 6X22 wearable configuration. Results are 
obtained by averaging output from 100 Monte Carlo simulations. 


6. Effect of Node Failures 


Node failures can occur when a sensor’s battery runs out or the sensor may be 


asleep when the source is detected. The scenario representing twelve soldiers with two 


a 


sensors each (6X22) is studied. Each soldier is stationed 10 m apart, and the source is 
100 m from the origin. Figure 31 shows the RMS error as a function of the number of 


sensor node failures. 


The RMS error ysteadily increases as the number of sensor node failures 


increases from two to twelve. For example, the RMS error for the Rs value increases 
gradually from 10 m to 20 m when the number of sensor node failures increases from 
zero to ten. This demonstrates the advantage of having distributed wireless sensor 
networks where sensor node failures lead to a graceful degradation of performance and 


not a sudden system failure. Note that drastic increase in the RMS error y at the SI stage 


is observed only when there are twelve or more sensor node failures (i.e., 50 % sensor 
node failure). At such a high failure rate, the SI method used in the first stage of the 
hybrid SI/ML method gives a large RMS error of 50 m for the Rs value, but the ML 
method in the second stage is able to tolerate the errors and produced RMS errors for the 


Rs value of less than 10 m. 
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Figure 31. RMS error vy as the function of the number of sensor node failures, 


representing a 6X22 wearable configuration. Results are obtained by averaging 
output from 100 Monte Carlo simulations. 
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Te Effect of Location Perturbations 

In field deployment, it is expected that there will be some perturbation (or slight 
deviation) in the sensor locations as perfectly straight line troop formation or exact 
symmetry is usually not possible. This perturbation effect on location estimation 
performance is studied by introducing some degree of randomness to the sensor 
positions. The standard deviation of the random distance between soldiers is varied from 
0.01 to 1; the standard deviation of the two sensor locations on a soldier is set to 0.01. 


Figure 32 shows the RMS error yfor different values of standard deviations of location 


perturbations. 


It is seen that RMS error vis similar for low and high standard deviations of 
location perturbation. For example, the RMS error for the Rs value is between 7 m to 9 
m for the ML method at the second stage of the hybrid SI/ML algorithm for the various 
levels of location perturbation standard deviations. This shows that location perturbation 
has limited effect on the estimation accuracy. This clearly illustrates the flexibility of the 


network configuration and demonstrates the suitability for potential field deployment. 
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Figure 32. RMS error for different values of standard deviations of location 


perturbations, representing a 6X2x2 wearable configuration 
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Cc; WIRELESS NETWORK SIMULATION 

The wireless network simulation is carried out using NS-2 simulator, a discrete 
event simulator suitable for networking research. The wireless extension and the 802.15.4 
extension [19] are used to conduct the studies on scenarios using an IEEE 802.11-based 
wireless local area network (WLAN) and an IEEE 802.15.4-based wireless personal area 
network (WPAN). Figure 33 shows the protocol architecture of the NS-2 wireless 


extension to support WPAN simulation. 






Wireless Scenario 
Definition 
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Upper Layers 


Routing 











Figure 33. The protocol architecture of the IEEE 802.15.4 WPAN extension in NS-2 
(from [8]). 


A scenario of a wireless sensor network consisting of 26 sensor nodes situated 
over a 50 mx50 m area is simulated using NS-2. Figure 34 shows the layout of the 
sensors defined for the scenario used for the simulation. The traffic scenario of twelve 
sensor nodes (#1 to #12) sending Constant Bit Rate (CBR) traffic to a central node (#19) 


at the same time is simulated. The Adhoc On-demand Distance Vector (AODV) routing 
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protocol [20] is run on top of the 802.15.4 MAC and PHY layers to study the 
performance of WPAN networks. Similarly, the AODV routing protocol is run on top of 
802.11 MAC and PHY layer to study the performance of WLAN networks. 
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Figure 34. NS-2 WPAN and WLAN Simulation Scenario: 26 sensor nodes in a 50 
x50 m’ area. The Time Interval Length (TIL) used in the simulation is equal to 
one second. 


i, Sending/Receiving Throughput 

Figure 35 shows the sending and receiving throughput for WPAN and WLAN 
networks. The Time Interval Length (TIL) used in the simulation is equal to one second. 
The average sending throughput refers to the average throughput of nodes 1 to 12 that are 
sending CBR traffic to central node 19; the average receiving throughput refers to the 


average throughout received by node 19. 


For WPAN networks, the average sending throughput is approximately 100 kbps 
and the receiving throughput is approximately 20 kbps. It is observed that the WLAN 


aT 


network has better throughput than WPAN network. For WLAN networks, the average 
sending throughput is approximately 300 kbps and the receiving throughput is 
approximately 200 kbps. 
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Figure 35. Sending and Receiving Throughput of: (a) WPAN and (b) WLAN. The 


Time Interval Length (TIL) used in the simulation is equal to one second. NS-2 
simulation of a 26-node network operating over an area of 50 x 50 m. 
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2. End-to-End Delay 

Figure 36 shows the histograms and the summary statistics for the end-to-end 
delay for the WPAN and WLAN networks. For WPAN networks, approximately 90 % of 
the traffic has end-to-end delay of less than 10 ms. The average WPAN end-to-end delay 


is 6.4 ms. The minimum WPAN delay is 0.35 ms and the maximum delay is 115 ms. 


It is observed that the end-to-end delay in the WLAN network is slightly higher 
than that in the WPAN network. For WLAN networks, approximately 97 % of the traffic 
has end-to-end delay of less than 50 ms. The average WLAN end-to-end delay is 17 ms. 
The minimum WLAN delay is 0.2 ms and the maximum delay is 564 ms. This 
observation of lower end-to-end delay in the WPAN networks compared to the WLAN 
networks is consistent with the results of the WPAN performance study performed by 


Zheng [19]. 
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Figure 36. End-to-End delay for: (a) WPAN and (b) WLAN. NS-2 simulation of a 


26-node network operating over an area of 50 x 50 m’. 


3s Packet Size 

Figure 37 shows the simulation summary statistics for WPAN and WLAN 
networks. The average WPAN packet size is 51.37 bytes with a minimum packet size of 
5 bytes and a maximum packet size of 111 bytes. The WLAN networks have a similar 
packet size as WPAN networks. The average WLAN packet size is 57.51 bytes with a 


minimum packet size of 28 bytes and a maximum packet size of 142 bytes. 
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Simulation information: Simulation information: 


Simulation length in seconds: 58.99960653 Simulation length in seconds: 58.86696361 
Number of nodes: 26 Number of nodes: 26 
Number of sending nodes: 26 Number of sending nodes: 26 
Number of receiving nodes: 23 Number of receiving nodes: 23 
Number of generated packets: 23206 Number of generated packets: 42689 
Number of sent packets: 21326 Number of sent packets: 42250 
Number of forwarded packets) 0 Number of forwarded packets: 0 
Number of dropped packets: 16527 Number of dropped packets: 7929 
Number of lost packets: 6866 Number of lost packets: 11339 
Minimal packet size: 5 Minimal packet size: 28 
Maximal packet size: 11 Maximal packet size: 142 
Average packet size: 51.3744 Average packet size: 57.5045 
Number of sent bytes: 1151295 Number of sent bytes: 2635686 
Number of forwarded bytes: Oo Number of forwarded bytes: 0 
Number of dropped bytes: 1059966 Number of dropped bytes: 478598 





Packets dropping nodes: fo 12345678910 Packets dropping nodes: fo Wie él ef EP 


(a) (b) 


Figure 37. Simulation summary statistics for: (a) WPAN and (b) WLAN networks. 
NS-2 simulation of a 26-node network operating over an area of 50 x 50 m’. 


4. Jitter 
Figure 38 shows the jitter plots of transmitted and received traffic for WPAN and 
WLAN networks. It is observed that the packet jitter is low for both WLAN and WPAN 


networks. 


For the WPAN, the maximum jitter of sending packets is 200 ms and the 
maximum jitter of receiving packets is 320 ms. The average jitter of sending and 
receiving packets is approximately 8 ms. For WLAN, the maximum jitter of sending 
packets is 200 ms and the maximum jitter of receiving packets is also 200 ms. The 


average jitter of sending and receiving packets is approximately 5 ms. 
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Figure 38. Jitter plot of: (a) all sent WPAN packets, (b) all received WPAN packets 
(c) all sent WLAN packets, (d) all received WLAN packets. NS-2 simulation of a 
26-node network operating over an area of 50 x 50 m. 


The simulation results show that the wireless sensor network has a throughput 
between 20 kbps and 100 kbps, a low average delay of less than 10 ms and a low average 


jitter in the range of 5 ms to 10 ms. Assuming that the payload sent by each sensor to a 
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server for processing is 1 kByte, as has been implemented by Duckworth et al. [6], the 
required transmission time is approximately between 0.08 s and 0.4 s. After considering 
the processing time and other delays, the total delay can be approximated to be less than 
one second. This indicates that the wireless sensor network is suitable for implementation 


in fast-response counter-sniper applications. 


This chapter has presented the simulation results of the hybrid SI/ML estimation 
method used in the two-step localization process for counter-sniper applications. 
Simulation results indicate that the proposed two-step localization method is suitable for 
counter-sniper applications. For example, twelve soldiers with two sensors each (6X2 x2) 
can detect a sniper 100 m away with an RMS error of 3 m, 8 m and 15 m for the Xs, Ys 
and Zs values, respectively. Better accuracy can be obtained by increasing the number of 
sensors or by increasing the inter-sensor spacing. The NS-2 based simulation results 
indicate that both WPAN and WLAN networks are suitable for implementation in fast- 
response counter-sniper applications due to their high throughput, low delay and low 
jitter. The next chapter will conclude the thesis by summarizing the significant results and 


proposing possible future work. 
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V. CONCLUSION 


This thesis focused on source localization using wireless sensor networks. 
Wireless sensor nodes worn on soldiers or installed on vehicles can be deployed to form a 
distributed sensor array to locate the source of sniper fire. A two-step source localization 
process is proposed for counter-sniper applications: the time difference of arrival 
(TDOA) values are first determined using the generalized cross correlation (GCC) 
method and then the TDOA values are used by a hybrid spherical interpolation/maximum 
likelihood (SI/ML) estimation method to determine the shooter location. A simulation 
model has been developed in MATLAB to study the performance of the hybrid SI/ML 
estimation method. A wireless sensor network is simulated in NS-2 to study the network 


throughput, delay and jitter of the IEEE WPAN and WLAN networks. 


A. SIGNIFICANT RESULTS 

The effects of number of sensor nodes, inter-sensor spacing, topology, noise, 
source-sensor distance, sensor node failure, and location perturbation on the location 
estimation performance of the hybrid SI/ML method have been studied through 
MATLAB simulation. 


Simulation results show that the error variance can be reduced by increasing the 
number of sensor nodes or the inter-sensor spacing. To increase the estimation accuracy, 


the number of sensors as well as the inter-sensor spacing should be increased. 


For the case of wearable sensors, the constraint of small inter-sensor spacing on 
the human body degrades the estimation accuracy. On the other hand, vehicular 
configuration providing larger vertical and horizontal inter-sensor spacing can help 
improve the estimation accuracy. We see that the configuration of a group of twelve 
soldiers or vehicles, with two wireless sensors per soldier or vehicle, is sufficient to 
accurately detect a sniper situated 100 m away with an RMS error of 3 m, 8 m and 15 m 
for the Xs, Ys and Zs values, respectively. By using a longer range wireless network such 


as WLAN, the inter-sensor spacing can be increased and provides better estimate 


65 


accuracy to detect a sniper located farther. Simulation results also indicate that the sensor 


topology should be well represented in all three dimensions to obtain desirable resolution. 


The source estimation accuracy is not adversely affected when encountering node 
failures or location perturbations. Results show that sensor node failures lead to a gradual 
degradation of accuracy and not a sudden total system failure. This demonstrates the 
advantage of using distributed wireless sensor networks to provide greater system 
robustness. It is observed that location perturbations do not adversely affect the location 
estimation. This illustrates the flexibility of the network configuration and demonstrates 


the suitability for potential field deployment. 


Network simulation results based on the NS-2 package show that the wireless 
sensor network can support a throughput of between 20 kbps and 100 kbps, a low average 
delay of less than 10 ms, and a low average jitter in the range of 5 ms to 10 ms. These 
are suitable for fast response network-based sniper detection applications that require low 


delay and fast data exchange. 


B. FUTURE WORK 

This thesis dealt with a single source scenario. In a future effort, the work 
reported here can be extended to multiple sources. Since the SI method is only applicable 
to a single source case, a different estimation method should be considered. Other source 


types, such as RF and seismic signals, can also be included in the simulation model. 


We have seen that the SI method provided good initial values for the ML method. 
However, the SI method can introduce bias at low noise levels and affect the estimation 
accuracy of the ML method. Alternative estimation methods, such as linear interpolation 
(LD) [21] and cluster-based schemes [16], can be studied to replace the SI method in the 
two-step process. We also observed that the ML method failed to converge when the SI 
method produced noisy location estimates. Further study is required to analyze this 


behavior. 


To improve the estimation accuracy of a nearby source, beamforming techniques 
can be applied as they generally provide a more accurate and robust performance in the 


near-field [17][22][23]. As discussed in this thesis, the proposed two-step process is 
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suitable for far-field applications. A combination of beamforming and the two-step 


process can be studied in a future effort to cover both near-field and far-field sources. 


In this work, the two-step source localization process was simulated in MATLAB. 
This work can be extended by implementing the scheme in hardware using Crossbow 
motes and conducting field experiments. The measured results can be used to validate the 


simulated results provided in this thesis. 
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APPENDIX: MATLAB SOURCE CODE 


This appendix lists the MATLAB source code used to study the performance of 
the Spherical Interpolation (SI) and Maximum Likelihood (ML) source localization 


estimation method. 


% Calculate Source Location based on Range Difference of Arrival (RDOA) 
% (1) Least Square Estimation - Spherical Interpolation (SI) Method 

% (2) Maximum Likelihood (ML) Method 

% 

% N sensors, 1 source 

% Using sensor | as reference i.e (x1=0, y1=0, z1=0) 

% 


nRun=100; % number of Monte Carlo runs 
% uncomment one of them 

% bML=0; % turn off ML calculation 
bML=1; % turn on ML calculation 


% uncomment one of them 


perturb=0; % turn off location perturbation 
% perturb=1; % turn on location perturbation 
ee ae eae a eee eee eee ee ee ee ae ee ee 


% Actual source location (m) in Cartesian coordinates x, y and z 
% Note: For simplicity, we only varies y for our simulation 

OO Se Se ae ee 
xs_src_actual=[0]; 


% Varies the Y position (Choose 1) 


ys_src_actual=[100]; %100 m 
zs_src_actual=[0]; 


Rs_actual=sqrt(xs_src_actual.’2 + ys_src_actual.2 + zs_src_actual.’2); 
% calculate corresponding range Rs 

bearing_actual=[xs_src_actual; ys_src_actual; xs_src_actual]/Rs_actual; 
% calculate corresponding bearing 
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% Actual sensor location (m) in Cartesian coordinates x, y and z 

% Note: For simplicity, we only use integers and then multiply with 
% a scaling factor to produce the actual coordinates. 

% e.g. [5 10 15] =[1 23] *5; 

oe ee Ona Ree Ae Ee ee eae ee eee oe eee 


scale_dist = 10 % 10m 
% (Choose 1 of the following sensor configuration for study) 


% 12x2 sensors arranged 2 rows 

% xi=[(001122334455 001122334455 ].*scale_dist; 

% yi=[0 10101010101 01010101010 1).*scale_dist; 

% 2i=[(000000000000 11111111111 1).*1.0; 

% Soldier configuration (each with 2 sensors). z= 1m apart vertically 


temp=size(x1); 
nSen=temp(1,2); % number of sensor (>4) 


noisestd=1; 
if (perturb==1) 
randn('state',0); 
tmp1=randn(3, nSen); 
for i=1:2:nSen 
xi(i)= xi(i1) + noisestd*tmp1(1,1); 
xi(i+1)= xi(i+1)+(noisestd+0.01)*tmp1(1,1); % less variance on the body 
yi(i)= yi(i) + noisestd*tmp1 (2,1); 
yi(it+1)= yi(@i+1)+(noisestd+0.01)*tmp1(2,1); % less variance on the body 
end 
zi=zi+0.01*tmp1(3,:); % less variance on the body 
end 


% RD noise (Choose 1) 

ae Tet UA OES TALE ERATE 

Noise_Factor=0.001; % noise std = Std_Norm * (source distance). % 
we expect bigger noise variance for larger distance. 
Noise_Var=(Noise_Factor*Rs_actual)2; 
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% Random Process 

% AWGN 
randn(‘state',Q); 
noise = sqrt(Noise_Var)*randn(nRun, 1); 
Ynoise_mean = mean(noise, 2); % average along row 


for k=1:nRun, % Monte Carlo Simulation 


Ai=[x1 yi zi |; 
Di= sqrt ((xi-xs_src_actual).“2 + (yi-ys_src_actual).42 + (zi-zs_src_actual).‘2); 
Ri= sqrt ((xi).42 + (yi).42 + (zi).42); 


locSen=[x1' yi' zi']; 


% using N sensors 
for i=1:nSen-1 
%d21=Di(2)-Di(1); 
%d31=Di(3)-Di(1);.. 
%dn1=Di(n)-Di(1); 
%d=[d21;d31;...;dn1]; 
d(i, 1)=Di(i+1)-Di(1)+noise(k); %add noise to RD estimates 


% delta2=Ri(2)*2-d(1)%2; 

% delta3=Ri(3)*2-d(1)2.... 

% deltan=Ri(n)*2-d(1)42; 

% delta=[delta2;delta3;...deltan]; 
delta(i, 1)=Ri(i+1)*2-d(1)2; 


% s2= [xi(2) yi(2) zi(2)]; 

% s3= [x1(3) yi(3) zi(3)];... 

% sn= [xi(n) yi(n) zi(n)]; 

% s=[s2;s3;...sn]; 
s(i,:)=[x10i+1) yiG+1) ziG+1)]; 


end 


% define weight (positive definite, i.e., diagonal positive and symmetrical) 
w=eye(nSen-1); % set to identity matrix for unweighted case 


Sw=(s'*w*s)\(-1)*s'*w; 
Ps=s*Sw; 
Ps_ortho=eye(nSen-1)-Ps; 


Rs_SI_cal=0.5*(d'*Ps_ortho*w*Ps_ortho*delta)/(d'*Ps_ortho*w*Ps_ortho*d); 


% Calculate Xs for SI method 

Xs_row_SI = 0.5*Sw*(delta-2*Rs_SI_cal*d); 
Xs_SI(k,:)=Xs_row_SI' ; 

Rs_SI(k,:)=sqrt(Xs_SI(k,1)42 + Xs_SI(k,2)42 + Xs_SI(k,3)%2); 


bearing_SI(k,:)=Xs_SI(k,:)/Rs_SI(k,:); 
% error_row= delta - 2*Rs_ SI*d - 2*s*Xs_row_SI; %error 
% error(k,:)=error_row'; 


% % Maximum Likelihood Method 
% % Objective function is contained in mlobjfun.m 


if (bML==1) 
x0 = Xs_SI(k,:);  % As value obtained from SI as starting guess 
% x0 = [0 ys_src_actual 0]; % Starting guess 
options = optimset('LargeScale’,'off’); 


% % LevenbergMarquardt 
options=optimset(options,'LevenbergMarquardt','on'); % LM 
% options=optimset(options, 'LevenbergMarquardt','off'); % Gauss Newton 
[x,resnorm,residual,exitflag,output]= 
Isqnonlin( @ mlobjfun,x0,[],[],options,locSen,Noise_Var,d); 

Xs_ML(k,:)=x; 
Rs_ML(k,:)=sqrt(Xs_ML(k, 1)42+Xs_ML(k,2)42+Xs_ML(k,3)%2); 
bearing_ML(k,:)=Xs_ML(k,:)/Rs_ML(k,:); 

end 


% Calculate bias (i.e., errors) for source location, range and bearing 
% estimates 


% SI 

bias_Xs_SI(k,1)=Xs_SI(k,1)-xs_src_actual; 
bias_Xs_SI(k,2)=Xs_SI(k,2)-ys_src_actual; 
bias_Xs_SI(k,3)=Xs_SI(k,3)-zs_sre_actual; 


% ML 
if (bML==1) 
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bias_Xs_ML(k,1)=Xs_ML(k,1)-xs_src_actual; 
bias_Xs_ML(k,2)=Xs_ML/(k,2)-ys_src_actual; 
bias_Xs_ML(k,3)=Xs_ML(k,3)-zs_src_actual; 
end 
end 


bias_Rs_ SI=Rs_SI-Rs_actual; 
bias_bearing_SI = 180/pi*acos(bearing_SI*bearing_actual); 


if (bML==1) 

bias_Rs ML=Rs ML-Rs_actual; 

bias_bearing_ML = 180/pi*acos(bearing_ML*bearing_actual); 
end 


meanxs_SI=mean(bias_Xs_SI(:,1)); 
meanys_SI=mean(bias_Xs_SI(:,2)); 
meanzs_SI=mean(bias_Xs_SI(:,3)); 
meanrs_SI=mean(bias_Rs_ SI); 
meanbear_SI=mean(bias_bearing_ SI); 


vect_mean_SI=[meanxs_SI;meanys_SI;meanzs_SI;meanrs_SI;meanbear_SI]; 


%ML 

if (bML==1) 

meanxs_ML=mean(bias_Xs_ML(:,1)); 

meanys_ML=mean(bias_Xs_ML(:,2)); 

meanzs_ML=mean(bias_Xs_ML(:,3)); 

meanrs_ML=mean(bias_Rs_ ML); 

meanbear_ML=mean(bias_bearing_ML); 

vect_mean_ML=[meanxs_ML;meanys_ML;meanzs_ML;meanrs_ML; 
meanbear_ ML]; 

end 


% Calculate Variance = E[(a - mean)*2] 

a ee eS Cee ee a AL EO 
varxs_SI=var(bias_Xs_SI(-:,1)); 
varys_SI=var(bias_Xs_SI(:,2)); 
varzs_SI=var(bias_Xs_SI(:,3)); 
varrs_SI=var(bias_Rs_ SI); 
varbear_SI=var(bias_bearing_ SI); 


vect_var_SI=[varxs_SI;varys_SI;varzs_SI;varrs_SI;varbear_S1]; 


%ML 
if (bML==1) 
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varxs_ML=var(bias_Xs_ML(:,1)); 
varys_ML=var(bias_Xs_ML(:,2)); 
varzs_ML=var(bias_Xs_ML(:,3)); 
varrs_ML=var(bias_Rs_ ML); 
varbear_ML=var(bias_bearing_ML); 


vect_var_ML=[varxs_ML;varys_ML;varzs_ML;varrs_ML;varbear_ML]; 
end 


% Calculate second moment (RMS)= sqrt {E[a*2]} = sqrt {mean’2 + variance } 
Oe cts oe ee at Rhea 
rmsxs_SI=sqrt(mean(bias_Xs_SI(:,1)).42+varxs_SI); 
rmsys_SI=sqrt(mean(bias_Xs_SI(:,2)).42+varys_SD; 
rmszs_SI=sqrt(mean(bias_Xs_SI(:,3)).42+varzs_SD; 
rmsrs_SI=sqrt(mean(bias_Rs_SI).42+varrs_SI); 
rmsbear_SI=sqrt(mean(bias_bearing_SI).’2+varbear_SD; 


vect_rms_SI=[rmsxs_SI;rmsys_SI;rmszs_SI;rmsrs_SI;rmsbear_SI]; 


%ML 

if (bML==1) 
rmsxs_ML=sqrt(mean(bias_Xs_ML/(:,1)).42+varxs_ML); 
rmsys_ML=sqrt(mean(bias_Xs_ML/(:,2)).42+varys_ML); 
rmszs_ML=sqrt(mean(bias_Xs_ML(:,3)).42+varzs_ML); 
rmsrs_ML=sqrt(mean(bias_Rs_ML).42+varrs_ML); 
rmsbear_ML=sqrt(mean(bias_bearing_ML).42+varbear_ML); 


vect_rms_ML=[rmsxs_ML;rmsys_ML;rmszs_ML;rmsrs_ML;rmsbear_ML]; 
end 


% Calculate Cramer Rao Bound 
% 


cov_mat=Noise_Var.*(0.5*ones(length(d))+0.5*eye(length(d))); 


for i=1:length(d) 
al=[xs_src_actual-locSen(i+1,1) ys_src_actual-locSen(i+1,2) 
zs_src_actual-locSen(i+1,3)]; 


a2=sqrt((xs_src_actual-locSen(i+1,1))*2+(ys_src_actual-locSen(i+1,2))42 
+(zs_src_actual-locSen(i+1,3))2); 

bl=[xs_src_actual-locSen(1,1) ys_src_actual-locSen(1,2) 
zs_src_actual-locSen(1,3)]; 

b2=sqrt((xs_src_actual-locSen(1,1))42+(ys_src_actual-locSen(1,2))42 
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+(zs_src_actual-locSen(1,3))*2); 
jacobian(i,:)= (al/a2)-(b1/b2); 
end 


fisher=jacobian'*inv(cov_mat)*jacobian; 


crlb= trace(fisher’-1); % compare with MSE of Rs 


% hfig1=figure; 

if (bML==1) 

plot(xi, yi,'kv', xs_src_actual, ys_src_actual, 'k*', Xs_SI(:,1), Xs_SI(,2), 
‘ko'".Xs_ML(:,1), Xs_ML(.,2), 'kd'); % plot both SI and ML 

else 

plot(xi, yi,'kv', xs_src_actual, ys_src_actual, 'k*', Xs_SI(:,1), Xs_SI(-,2), 'ko'); 
% plot just SI only 

end 


title('Sensor and Source Location’); 

strl=sprintf(‘[Xs, Ys, Zs, Rs, Bearing], Noise Std = %s*Rs',Noise_Factor); 
str2=sprintf(‘SI Method’); 

str3=sprintf((RMS = [%s, %s, %s, %s, Yos]', rmsxs_SI, rmsys_SI, rmszs_SI, 
rmsrs_SI, rmsbear_SI); 

str4=sprintf(‘Mean = [%s, %s, %s, %s, %s]', meanxs_SI, meanys_SI, meanzs_SI, 
meanrs_SI, meanbear_SJ); 

strS=sprintf('Variance = [%s, %s, %os, %s, %os]', varxs_SI, varys_SI, varzs_SI, 
varrs_SI, varbear_SI); 

if (bML==1) 

str6=sprintf(ML Method’); 

str7=sprintf((RMS = [%s, %s, %s, %s, %os]', rmsxs_ML, rmsys_ML, rmszs_ML, 
rmsrs_ML, rmsbear_ML); 

str8=sprintf(‘Mean = [%s, %s, %s, %s, Yos]', meanxs_ML, meanys_ML, 
meanzs_ML, meanrs_ML, meanbear_ML); 

str9=sprintf('Variance = [%s, %s, Ys, Ys, %os]', varxs_ML, varys_ML, 
varzs_ML, varrs_ML, varbear_ML); 

str=sprintf('%s \n%s \n%s \n%s \n%s \n%s \n%s \n%s \n%s', str1, str2, str3, str4, 
str5, str6, str7, str8, str9); 

legend(‘sensor location’, ‘actual source location ', ‘calculated source location 
(SD','calculated source location (ML)'); 

else 

str=sprintf('%s \n%s \n%s \n%s \n%s \n%s \n%s \n%s \n%s', str1, str2, str3, str4, 
str5); 

legend(‘sensor location’, ‘actual source location’, 'calculated source location (SI)'); 
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end 


text(scale_dist/10,ys_src_actual,str); 
xlabel('‘Distance (metres) in X direction’); 
ylabel('Distance (metres) in Y direction’); 


% generate results output files 
fid = fopen(‘results.txt’,'w’); 


for k=1:nRun, 

fprintf(fid,'%e\t%e\t%e\t%e\t%e\n' bias_Xs_SI(k,1),bias_Xs_SI(k,2), 

bias_Xs_SI(k,3), bias_Rs_SI(k), bias_bearing_ SI(k)); 

end 
fprintf(fid,'\n%e\t %e\t Ze\t %e\t %e\n', meanxs_SI, meanys_SI, meanzs_SI, 
meanrs_SI, meanbear_SI); 
fprintf(fid,'%e\t %e\t %e\t Z%e\t %e\n', varxs_SI, varys_SI, varzs_SI, varrs_SI, 
varbear_SI); 
fprintf(fid,'%e\t %e\t %e\t %e\t %e\n', rmsxs_SI, rmsys_SI, rmszs_SI, rmsrs_SI, 
rmsbear_SI); 


fclose(fid); 
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